#include <iostream>
#include<sys/time.h>
#include<arm_neon.h>
#define N 1000
#define C 1
using namespace std;
float A[N][N];


void init0(float A[N][N]) {
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            A[i][j] = 0;
        }
        A[i][i] = 1.0;
        for(int j = i + 1; j < N; j++) {
            A[i][j] = rand();
        }
    }
    for(int k = 0; k < N; k++) {
        for(int i = k + 1; i < N; i++) {
            for(int j = 0; j < N; j++) {
                A[i][j] += A[k][j];
            }
        }
    }
}//initialize

void normal(float A[N][N]) {
    struct timeval h, t;
    float time = 0.0;
    for(int r = 0; r < C; r++) {
        init0(A);
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                for(int j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME1: "<<time/C<<" ms"<<endl;
}



//Neon全部并行（不对齐）
void neon(float A[N][N]) {
    struct timeval h, t;
    float time = 0.0;
    for(int r = 0; r < C; r++) {
        init0(A);
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            float* akk = new float[4];
            *akk = *(akk+1) = *(akk+2) = *(akk+3) = A[k][k];
            float32x4_t vt = vld1q_f32(akk);
            delete[] akk;
            for(int j = k + 1; j + 3 < N; j += 4) {
                float32x4_t va = vld1q_f32(&A[k][j]);
                va = vdivq_f32(va, vt);
                vst1q_f32(&A[k][j], va);
            }
            for(int j = N - N%4; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = A[i][k];
                float32x4_t vaik = vld1q_f32(aik);
                delete[] aik;
                for(int j = k + 1; j + 3 < N; j += 4) {
                    float32x4_t vaij = vld1q_f32(&A[i][j]);
                    float32x4_t vakj = vld1q_f32(&A[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);
                    vaij = vsubq_f32(vaij, vx);
                    vst1q_f32(&A[i][j], vaij);
                }
                for(int j = N - N%4; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME2: "<<time/C<<" ms"<<endl;
}
//Neon第一处并行
void neon1(float A[N][N]) {
    struct timeval h, t;
    float time = 0.0;

    for(int r = 0; r < C; r++) {
        init0(A);
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            float* akk = new float[4];
            *akk = *(akk+1) = *(akk+2) = *(akk+3) = A[k][k];
            float32x4_t vt = vld1q_f32(akk);
            delete[] akk;
            for(int j = k + 1; j + 3 < N; j += 4) {
                float32x4_t va = vld1q_f32(&A[k][j]);
                va = vdivq_f32(va, vt);
                vst1q_f32(&A[k][j], va);
            }
            for(int j = N - N%4; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                for(int j = k + 1; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME3: "<<time/C<<" ms"<<endl;
}
//Neon第二处并行
void neon2(float A[N][N]) {
    struct timeval h, t;
    float time = 0.0;
    for(int r = 0; r < C; r++) {
        init0(A);
        gettimeofday(&h, NULL);
        for(int k = 0; k < N; k++) {
            for(int j = k + 1; j < N; j++) {
                A[k][j] = A[k][j]/A[k][k];
            }
            A[k][k] = 1.0;
            for(int i = k + 1; i < N; i++) {
                float* aik = new float[4];
                *aik = *(aik+1) = *(aik+2) = *(aik+3) = A[i][k];
                float32x4_t vaik = vld1q_f32(aik);
                delete[] aik;
                for(int j = k + 1; j + 3 < N; j += 4) {
                    float32x4_t vaij = vld1q_f32(&A[i][j]);
                    float32x4_t vakj = vld1q_f32(&A[k][j]);
                    float32x4_t vx = vmulq_f32(vaik, vakj);
                    vaij = vsubq_f32(vaij, vx);
                    vst1q_f32(&A[i][j], vaij);
                }
                for(int j = N - N%4; j < N; j++) {
                    A[i][j] = A[i][j] - A[i][k]*A[k][j];
                }
                A[i][k] = 0;
            }
        }
        gettimeofday(&t, NULL);
        time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    }
    cout<<"TIME4: "<<time/C<<" ms"<<endl;
}


int main()
{

        normal(A);

        neon(A);

        neon1(A);

        neon2(A);


    cout<<endl;
    return 0;
}